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We analyze the phase diagram of a model of hard spheres of chemical radius one, which is defined 
over a generalized Bethe lattice containing short loops. We find a liquid, two different crystalline, a 
glassy and an unusual crystalline glassy phase. Special attention is also paid to the close-packing limit 
in the glassy phase. All analytical results are cross-checked by numerical Monte-Carlo simulations. 

PACS numbers: 64.70.Pf,64.60.Cn,75.10.Nr 

I. INTRODUCTION 

Even after many years of vivid interest, the structural glass transition is still an open and alive topic of research 
Q IE ^ 1^' Being defined by a drastic slowing down of the equilibration time of a liquid, it is introduced as a 
dynamical, non-equilibrium phenomenon. Up to now it is, however, one of the crucial questions if this slowing down 
is accompanied by an equilibrium glass transition or not. 

Recently, lattice-gas models have played an increasing role in this discussion Q, S 13 ■ The central idea of these 
models is to incorporate static geometrical constraints via hard-sphere interactions which restrict the possible particle 
packings and introduce some kind of geometrical frustration. Defined on (generalized) Bethe lattices, these models 
are characterized by the existence of a dynamical glass transition, followed by a static one. 

These models are contrasted by the so-called kinetically constrained models , in which not all particle moves are 
permitted, but which are characterized by a trivial thermodynamic equilibrium behavior .Defined on a Bethe lattice, 
they show, however, a very similar dynamical behavior to the above mentioned models (T(]| . 

The model analyzed in this paper was introduced in , and it belongs to the class of geometrically constrained 
models. Its equilibrium behavior will be analyzed in great detail using the cavity method [ill . The model can 
be defined over any lattice or graph. Sites are either occupied by particles, or they are empty. The particles are, 
however, too large to allow any two neighboring sites to be occupied simultaneously. In the presence of short loops in 
the lattice, this hard-core exclusion is sufficient to create a very rich phase diagram, including a liquid, two different 
crystalline, a glassy and an unusual, crystalline glassy phase. Whereas this paper concentrates completely on the 
static properties of the model, a subsequent publication jl3j will consider the dynamics. 

The paper is organized as follows: In the next section, the model is defined. Sec. Ill is dedicated to a description 
of an iterative method for the calculation of the partition function, or more precisely of effective local fields. Sec. IV 
solves this iterative method for the liquid-crystal transition, whereas Sees. V and VI are dedicated to the glass phase. 
Conclusion and outlook are given in the last section. 

II. THE MODEL 

The model was already introduced and studied in 0. In its most general formulation, it is defined on a graph G 
having vertices V = {1, ■■■,N} and undirected edges {i,j} € E connecting pairs of vertices. Vertices can either be 
empty, Ui = 0, or they are occupied by a particle, Ui = 1. These particles interact via a hard core of chemical radius 
\ one, i.e. neighboring vertices cannot be occupied simultaneously by two particles. Formulated in a more formal way, 
all edges {i,j} G E fulfill the constraint ntUj = 0, i.e. one of the end vertices has to be empty. The model can thus 
be characterized by its grand-canonical partition function 

S(^)= ^ e'^Ef..". Jl (l-n,n,), (1) 

ni,...,njve{0,l} {»,j}e£ 



or its grand-canonical potential 



n = --InS. (2) 
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FIG. 1: Part of a Bethe lattice with k = 2. Black vertices mark occupied sites, gray ones cannot be occupied due to volume 
exclusion. The lower right branch shows the densest local packing. 



Here we have introduced the chemical potential /i which is coupled to the total particle number and can thus 

be used to regulate the particle density 



Note that the product over all edges in Eq. serves as an indicator function for allowed particle configurations: 
Whenever there is at least one pair of occupied adjacent vertices, the product vanishes and thus the configuration 
does not contribute to the partition function. Note also that due to the hard-core interaction between particles the 
temperature does not play any role in the model, without loss of generality it is set to one. 

So far, the model is still defined on a general graph. In this work we will concentrate on (generalized) Bethe lattices. 
Bethe lattices can be defined in two different ways: 

• The first definition considers a Bethe lattice as an infinite regular tree, i.e. a cycle-free graph of fixed vertex 
degree (or coordination number) fc -|- 1. The hard-sphere model on this graph was first considered by Runnels 
0. 

• As a second possibility we can define a Bethe lattice as a random {k + l)-regular graph of N vertices, with 
A'' ^ oo in the thermodynamic limit. This version is frequently used in the context of glassy systems |l2l| . and 
allows in particular for the simulation of finite systems. 

Locally, both systems are equivalent due to the fixed vertex degree, the case k = 2 together with a possible particle 
packing is visualized in Fig. ^ The crucial difference results from the existence of cycles in random graphs: Even if 
these are of length O(lniV), i.e. their length diverges for large N, they are in general inconsistent with the crystalline 
structure which will be considered in Sec. IIVI A densest particle packing on a Bethe lattice defined as an infinite tree 
consists of an alternation of occupied and empty sites. In random regular graphs, many cycles of odd length exist, 
preventing this alternating configuration from being globally feasible. A way out of this dilemma will be discussed 
later on. 

Realistic, i.e. finite-dimensional, systems contain, however, many short cycles. To imitate this, we generalize Bethe 
lattices in the following way, cf. also Fig. [3 The graph is composed of cliques, i.e. fully connected subgraphs of 
p + I vertices each. In each vertex, fc -t- 1 of these cliques merge, such that the resulting graph has constant vertex 
degree p(k + 1). The global structure resembles, however, a Bethe lattice: We assume that besides the loops inside 
the cliques there are no other short cycles. Note that ordinary Bethe lattices are obtained in the special case p — 1, 
which is included in the following discussion of the general case. Also generalized Bethe lattices can be defined in the 
two ways discussed above, with similar consequences on the global feasibility of locally dense packings. 




(3) 
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FIG. 2: Part of a generalized Bethe-lattice with fe = 2, p = 3. The 21 cliques connecting the central vertex to it's nearest and 
second neighbors are shown. 




FIG. 3: Iteration of rooted trees: k = 2 branches of p = 3 rooted trees are merged to form a new tree with root 00. 



III. ITERATION OF THE PARTITION FUNCTION 

The direct calculation of the partition function is complicated. It can, however, be achieved via an iterative method. 
For doing so, we shall introduce the notation of a rooted subtree: Imagine one of the cliques containing a vertex i to be 
removed from the graph, then a rooted subtree consists of the connected component containing i. The vertex i itself 
is denoted as the root of the subtree, cf. Fig. |31 Note that the root has only kp neighbors, whereas all other vertices 
in the subtree have degree (fc + l)p. For each of these rooted trees we introduce two restricted partition functions: 
denotes the partition function with rii fixed to zero, whereas corresponds to an occupied root. 

Given kp rooted trees {j = 1, . . . , k, I = 1, . . . ,p), these can be composed according to Fig.|3|to form one new rooted 
tree with root j = 0, ^ = 0. The double index is introduced to lighten the notation in the following equations: The 
first index defines the (removed) clique containing the root, whereas k enumerates the vertices within this clique. The 
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restricted partition function of the new tree can now be easily inferred from those of the kp original ones: 

k / p p p \ 



n n-i'+E-^*'n-r (4i) 



= 1 ;=i m=l 



k 

"00 _ „M 



j=i / 

The expression in the brackets stands for the contribution of one branch: In case that the root 00 is empty, up to one 
of the neighbors in each branch can be occupied. If, on the other hand, the root is occupied, none of its neighbors 
is allowed to be occupied. The additional factor e'' in the second equation takes into account the chemical potential 
acting on the root particle. 

Introducing now the cavity fields 

(more precisely, ln(x^' + defines a local field coupled to n^/), we find 



= i(x",a;i2,...,a;'=f). (6) 

In order to understand the thermodynamic behavior of our model, we have to find appropriate solutions to these 
equations. In the next section, we shall discuss the case of the transition toward an ordered, crystalline phase, 
whereas Sec. is dedicated to glassy solutions. 

Before doing so, we have to clarify the physical meaning of x-'' and its connection to observables of the system. A 
simple observable is given by the local density p^" of vertex 00. This quantity can be obtained by composing rooted 
subtrees as well. In the original graph, each vertex has, however, (fc + l)p neighbors. Within one iteration step, this 
can be achieved by merging (fc + 1) instead of k branches. Denoting the corresponding partition functions by S^y^, 
we get 




^0" = TT f TT ^ Si' n sr j (7i) 

(7n) 



m— 1 



This results in 



^00 -^-oo 



(8) 



i.e. physical observables can be directly calculated from the cavity fields 



IV. CRYSTALLIZATION 



A. Liquid and crystalline solutions 

To solve Eq. @, we have to restrict the solution space by considering the limiting cases fi ±oo. For the moment 
we discuss only the first case for the definition of the generalized Bethe lattice, i.e. we exclude the existence of long 
loops which may have some influence on the global packing structure. 

The limiting case of an empty (or very dilute) system (/i — s- —oo) is characterized by a spatially homogeneous 
density describing a liquid phase. Consequently, the cavity fields are all equal to some x* given self-consistently by 



FIG. 4: Possible iteration steps leading to rooted trees with root in the 0-lattice (left), 1-lattice (right). Vertices belonging to 
the 0-lattice (1-lattice) are depicted by o (□). 



In the limit /i ^ oo we have to search for close-packings of the system. In these every clique carries exactly one 
particle, due to their regular structure these configurations are to be considered as being crystalline. The number of 
these configurations is exponential in the size of the graph (i.e. the number of sites), as can be seen easily: Initializing 
the configuration by putting one particle on an arbitrary vertex, all neighboring vertices have to be free. The {k + l)kp'^ 
second neighbors are partitioned into (fc -I- l)kp cliques oi p vertices each, and each of these cliques can carry exactly 
one particle - so there are possible configurations for selecting the particle positions in between the second 

neighbors. Iterating this argument, we obviously find an exponential number of close-packings. 

Selecting one close-packing, we can identify two sub-lattices: The first is formed by the occupied vertices (1-lattice), 
the second by the empty vertices (0-lattice). We introduce two cavity fields x*^*^^ and a;*^^^ for the sub-lattices. Note 
that a more general ansatz that works with p + I different cavity fields corresponding one to each site of a clique 
might lead to more complicated scenarios on (p + l)-partite graphs. It can be shown, however, that the cavity fields 
can only take two different values. We conclude that at least for the case p < 2 no extra solution appears. 

For the iteration of x*^"^ and x^^^ we have to distinguish two cases, represented in Fig.^ Either the root belongs to 
the 1-lattice, and has only neighbors from the 0-lattice, or it belongs to the 0-lattice, and has exactly one 1-neighbor 
and p — I 0-neighbors in each clique. The resulting equations thus read: 

=i(xW,x(°\...,a;(°),...) (9i) 

" V ' 

p-1 

=i(x("\...,x(°),...), (9ii) 

p 

where the presented block of variables has to be repeated fc times in the argument of x. The permutation symmetry 
of variables inside one block is already used here. The equations can be rewritten explicitly as 

= r (lOi) 

{i + xW + {p-i)x(o)y 

x(i) = r ■ (lOii) 

(l+px(o))' 



Using Eq. IHlwe thus find to sub-lattice densities 

p{0) 



(l + x(i) + (p-l)x(o)) 

/)(!) pA" 



fe+1 



(Hi) 

(llii) 



(l+px("))'+' ■ 

Eliminating the cavity fields using Eq. (|10|l . we find the equilibrium equations of state of our model: 

eA' = ^-^TTT 12i 

= ^ ' fc , 1 ■ 12ii 
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For fixed k and p, they imply functions p''^\^) and ^'^^•'(/i). Together with the grand-canonical potential which will 
be calculated in the following Scc. lIVBl they allow to determine the equilibrium behavior and, in particular, location 
and order of the crystallization transition. 

Let us therefore discuss the solutions of Eqs.El coming first back to the liquid solution. It can be understood as 
the special case of equal sub-lattice densities, p'"^ = p'^^ = p. We thus have to solve just one global equation of state, 



(l-(p+l)p) 

It has just one solution, growing monotonously from p — Q for /i — s- — oo to the maximally possible density p = 
for /i — > cx). 

Concentrating now on crystalline solutions with p*^°' ^ p^^\ we find that Eqs. (|12(l lead to the p-independent 
equilibrium condition 

p(^){l-pi^)f^pW{l-p(^))K (14) 

The function r(l — r)^ has a single maximum in r = -^p^^, thus for any p*-^) ^ ^^^-j- there is a solution for p^^^ 
being different from p(^\ only for p^^^ = both sub-lattice densities have to coincide. This defines a function 
p(o) = (/.(p(i)). Eqs. [121 can thus be reduced to the single equation 

pW(l-pW)'= 



(l-p(i)-p<?!)(p(i) ))'=+! ■ 

In order to find the spinodal point p^p, i.e. the minimal value of p where the last equation has a solution, we have to 
minimize the right-hand side of this equation with respect to p'^^\ 

For p''^^ ^ jjij- this leads to a condition on the sub-lattice densities exactly at the spinodal point, 

I- pW -pp(^) - pf){k + l){{kp~-l)p^^J -p+l) =0, (15) 

which is to be completed by the equilibrium condition p^^p = (j){p^^p). Taking as an example the case fc = 2, p = 2, 
there is exactly one solution to these equations: pip^ = i(6 — \/33) — 0.028 and pip = ^(9-|--\/33) ~ 0.819. The global 

density results in p^p = (ppip^ +Psp 1) — 0.292. From Eqs. (|12|l we infer the chemical potential p^p — 2.64 of the 

spinodal point. Plugging this into Eq. (|13f) in order to determine the liquid density at the same chemical potential, we 
find piiquid(psp) — 0.261. This value is substantially smaller than p^p, we thus conclude that the crystalline solution 
appears discontinuously at psp- Note also that, for p > p^p two different solutions for the sub-lattice densities are 
evolving, one with a monotonously increasing global density, the other one with an initially decreasing global density. 

For pip^ = jrpf, on the other hand, we know that pip^ — 0(pip^) = -j^-^ coincides with pip^ In this case, the 
crystalline solution appears continuously out of the liquid solution. Eqs. (|12|l are, however, consistent with this 
solution if and only if p = 1. The continuous onset of a crystalline solution can therefore be found only on usual 
Bethe-lattices. This statement does not include the scenario where an unstable crystalline solution becomes stable at 
p(o) = = _1_ which will be discussed later on. 



B. The grand-canonical potential 



In order to find out which of the solutions discussed above is the thermodynamically stable one, and if there 
are further metastable phases, we have to compare the corresponding values of the grand-canonical potential, or 
more precisely its density uj per clique. Due to the distinction of the two types of vertices, we can write this as 
Lo ~ puj'^^^ + u}^-^\ where 0;*^°/""^^ are the potentials per 0/1-site. Knowing the cavity fields x'"' and x'^^\ uj can be 
calculated using the following construction: 

• We start with p{k + 1) generalized Bethe-lattices. 

• We choose one clique in each lattice, and remove all its edges. 

• We obtain {p + 1) ■ p{k + 1) rooted trees, among which p ■ p{k + 1) have a 0-root, and p{k + 1) have a 1-root. 



7 



• We add p 0-vertices, and connect each one with fc + 1 of the 1-trees and {p — l){k + 1) of the 0-trees, such that 
we obtain p new generalized Bethe-lattices. 

• We add one 1-vertex, and connect it to the remaining p(k + 1) 0-trees, such that we get one more generahzed 
Bethe-lattice. 

• In this way, we obtain p + 1 lattices. 

We have to sum up the changes in the grand-canonical potential induced by following the construction: We have 
added p 0-vertices and one 1-vertex. To do so, we had to remove the edges of p{k + 1) cliques; denote the change 
in the potential induced by adding the edges of a clique by Af2iink. Further on we had to connect the newly added 
vertices and the corresponding rooted trees; the change in potential induced by one of these steps is denoted by Afig^j.^ 
resp. Ailgj?^^^ according to the added vertex. The total bilance thus reads: 



= + J'^ = -p{k + 1) Ar!n„k + P^n'^l + Anill . (16) 



Using definition Q of the grand-canonical potential, we calculate the above changes from the restricted partition 
functions of rooted trees: 



connected rooted trees 



e 



(sW-fS«)(S(°)-HSl°))^ 

V ' 

disconnected rooted trees 

1 -t-a;(i) +px''°'^ 



where we have divided both numerator and denominator by si^^ (Se'^'')P. Further on we get 

l-vertox occupied l-vertcx empty 



(17) 



free rooted trees 

_ e^-f(l+pa:W)fe+i 
~ (l-^a;(o))P('=+i) ' 

and 

(0, e^[Ei'\E^°^r-^]'^' + [Si^)(si°))^-i +Sl^)(si°))^-i + (p- l)si^)sl°)(s(°))^-^] 



(18) 



k+l 



[(sw + s«)(s(")-fsl°))p-i]^-+^ 

2^ + (l + .x(i) + (p-l)a;(0))''+' 



[{1 + xW){l + x(°))P^ 



k+l 



(19) 



This allows us to calculate the grand-canonical potential for the different crystalline and liquid solutions for the cavity 
fields, and thus to determine the phase diagram of the model by looking for the locally stable solution of minimal 
potential - provided the partitioning of the lattice in the two sublattices, i.e. the assumption of a crystalline structure, 
is consistent with the global structure of the lattice. 



C. The phase diagram 

Let us shortly summarize the case p = 1 where we have an ordinary Bethe lattice. We have found that at pcr = 
a crystalline solution appears continuously from the liquid one, i.e. we find a second-order crystallization transition. 
For fi < ficr, with p{fJ-cr) — Per, thc liquid solution is the only one, and thus globally stable. For fj, > pcr, this solution 
becomes locally unstable, and two stable crystalline states exist. They are connected to each other by exchanging the 
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two sub-lattices. For ^ oo, the system approaches one of its ground states which is given by a periodic alternation 
of occupied and empty sites. 

In the case of a generalized Bethe-lattice, i.e. for p > 1, the situation is slightly more involved. Below the spinodal 
point /isp only the liquid solution exists, it thus describes the globally stable solution. It exists, however, also for 
larger chemical potentials and, up to Hcr (as defined in the last paragraph), it remains at least locally stable under 
small perturbations. 

At fj,sp two other, crystalline solutions appear. The first one has always p'^^^ > p'"' , and the global density is larger 
than the liquid one. This solution is always locally stable, as can be shown by analyzing the relaxational dynamics of 
the system The second crystalline solution starts with the same sub-lattice densities at the spinodal point, but 
its global density first decreases. In addition it is initially locally unstable, hence unphysical. At /i = fj,cr, however, 
both sub-lattice densities become equal, p^-^^ = p^^^ — and exactly at this point it thus coincides with the liquid 
solution. Beyond pcr, the new solution becomes locally stable, and thus takes over the local stability from the liquid 
solution. Since both p(p)-cm:ves touch in this point with the same slope but different curvature, the corresponding 
transition there is of third order. The crystalline solution becomes also in some other aspect quite particular: The 
sub-lattice densities now fulfill p^^-* < p^^\ i.e. every clique has p vertices of higher, and one vertex of lower density! 
Note that this transition exists only if the global density can be reached. Every clique can be occupied by at most 
one particle, the density is thus restricted to values p < ^j^- A necessary and sufficient condition for the existence of 
the inverse crystalline solution is thus given by k > p. 

For all p > pLsp we thus find two locally stable solutions, and we have to consider the grand-canonical potentials in 
order to decide which solution is the globally stable one. For p > 1, the liquid solution remains the thermodynamically 
relevant one up to some /x^ > Psp, above which the first crystalline solution has the minimal w; a first-order freezing 
transition occurs. The second crystalline solution has never the globally minimal grand-canonical potential, but it 
crosses the one of the liquid solution at pcr, verifying thus the results of the local stability analysis. 

It has to be remarked that in the case of ordinary Bethe-lattices, i.e. p = 1, the equations are symmetric with 
respect to p^^^ and p^'^\ thus normal and inverse crystallization are indistinguishable. We can therefore consider this 
case as a degenerate one, where the crystalline and the inverse crystalline solutions appear in the same moment, and 
thus psp = Ps = Per holds. 

D. An example 

We want to illustrate this behavior in a specific example. The system of minimal k and p showing inverse crystal- 
lization has fc = 3 and p = 2. We therefore concentrate on these values. 

The phase diagram is shown in Fig.|Sl the grand-canonical potential per clique is shown in Fig.|Sl For small p, only 
the liquid phase exists. At psp ~ 1.386 a crystalline solution with sub- lattice densities pip^ ~ 0.034, p%^ ~ 0.637 and 
global density psp — 0.235 appears. Up to ps — 1.558, the liquid phase remains however the thermal equilibrium of 
our model. For all p above ps, the crystalline solution of higher density becomes the globally stable one, i.e. at ps 
the systems undergoes a first-order crystallization transition. The liquid solution stays, however, also locally stable 
beyond Ps, and in absence of sufficiently strong perturbations the system remains liquid also for higher p. Only 
at Per — 3.296 it loses its local stability in a third-order transition to an inverse crystalline state with p*^"^ > p^^^ 
Note that the density curves of the liquid and the second crystalline solution only touch at pcr, the inverse crystalline 
density is slightly larger than the liquid one, whereas the grand-canonical potentials cross each other with equal slopes. 

At this point we have to remark that the entropy s = —p{oJ -\- p) of the inverse crystalline solution vanishes at 
p ~ 5.58, and it becomes therefore unphysical at higher chemical potentials. Using a slight generalization of the 
IRSB formalism introduced in the next section, we find for p > 5.31 the existence of a solution that corresponds to 
a crystalline glass phase (cf. section IV A4ll which seems appropriate to solve the entropy crisis. Thus, the inverse 
crystallization is relevant only in an interval pcr < p < Peg with peg < 5.31. 

E. Comparison to numerical experiments 

To test our analytical results on the crystallization transitions, we have performed Monte-Carlo simulations. These 
are naturally performed on finite lattices, and a natural finite analog of the Bethe lattice defined as an infinite tree 
is provided by regular random graphs. As already discussed in Sec. m their large scale structure contains random 
loops which are not compatible with any division into 0- and 1-sub-lattices. This can be circumvented by imposing 
the consistency with the crystalline structure by generating, e.g., a random regular bipartite graph by grouping the 
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FIG. 5: Phase diagram for a generalized Bethe-lattice with fc = 3, p = 2. At a crystalline solution sets in discontinuously, 
and at /is it becomes globally stable. At ficr the liquid solution becomes locally unstable, and goes over into an inverse crystalline 
phase, which is relevant for /Xcr < n < Hcg with Hcg < 5.31. Locally stable solutions are denoted by full lines, locally unstable 
ones by dashed lines. 




FIG. 6: Grand-canonical potential per clique a; for fc = 3, p = 2. 



vertices into two subsets of sizes N and pN, and by introducing cliques only in between one vertex from the first and 
p vertices from the second subset. 

Note that also these graphs have an important difference to the generalized Bethe-lattices defined as infinite graphs 
which survives even in the thermodynamic limit: The first graphs allow only for exactly one partition into 0- and 
1-lattices, whereas the second definition allows for an exponential number of such partitions. This difference is 
meaningful for stability properties of the system. It can be shown in numerical simulations that even when going to 
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{p + l)-partite graphs, that aUow for p + 1 partitions, the inverse crystalhne phase becomes unstable. 

In the MC simulations we have included the following microscopic processes: Particles are mobile, i.e. they are 
allowed to move to neighboring sites whenever these are empty and do not have further occupied neighbors (hard- 
core constraint), and the system is coupled to a particle bath such that particles can be generated or deleted. The 
rates for the last two processes have to be related by detailed balance, i.e. they have ratio e'^, in order to guarantee 
equilibration. 

First, we have tested our predictions for the normal crystal phase on large graphs of = 6 • 10^ vertices. We 
have started the experiment in one densest packing, and then decreased the chemical potential with a rate of Sfi/St = 
—0.2/(100 MCs), and close to /isp with 6^/5t = —0.1/(200 MCs), starting with /i = 7. A further decrease in the 
decompactification rate did not alter the results, so we can expect thermal equilibrium. The left part of Fig. [71 shows, 
that the results are in perfect agreement with the analytical ones, in particular the discontinuous decrease of the 
density at the spinodal point is demonstrated with very high precision. 




0.35 - 



0.25 




0.15 



FIG. 7: Decompactification experiment for normal crystallization (left figure) and compactification experiment for inverse 
crystallization (right figure). Symbols which are rotated with respect to the legend correspond to smaller rates \5fi/St\, see 
main text. Analytical results are represented by full lines. 



The results for the inverse crystallization were tested by starting with an empty system which was compactified in 
a first step to /i = 1 where only the liquid phase and no crystalline phases exist. The system was compactified by 
increasing in discrete steps with two very small rates, 6^/6t — 0, 1/(10000 MCs) and 5fi/5t = 0, 1/(5000 MCs). 
Due to this slow dynamics we have also used a slightly smaller system of = 3 • 10^ vertices. 

The results are represented in the right of Figs.[7| The first observation is that, for small /i, numerical and analytical 
data coincide impressively: We find the liquid phase with equal sub-lattice densities for ^ < 3.3 and a transition to 
the inverse crystalline phase at ficr — 3.3, as predicted in the last section. For /i > ^cg — 4.1, the results depend 
on the compaction rate, underlining the glassy character of the corresponding phase. It is also observed that the 
global density is smaller than the inverse crystalline one, which mainly results from a lower sublattice density on the 
0-lattice. 

We mention that besides the transition from the liquid to the inverse crystalline phase which is shown in Fig. [7| 
(right), there is a certain probability for the system to undergo a transition to the normal crystalline phase with 
> p(o) "vvhen compacted beyond ficr- The dependence of this probability on the compaction rate and on the 
system size might be an interesting subject of further investigation. 



V. GLASSY BEHAVIOR 



In the last section, we have already seen that the simple crystalline ansatz of two cavity fields x*^"^ and x'^-* for 
solving the iterative Eq.El leads to inconsistencies for large fi. The entropy of the inverse crystalline solution, which 
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describes a metastable phase, becomes negative at finite /j,, and the solution is not physical any more. Frequently 
zero-entropy points can be understood as an indicator for the existence of a glass transition at some smaller value of 
the chemical potential, and thus a more general solution ansatz has to be considered. 

The situation is even more drastic if we consider random regular graphs as realizations of Bethe lattices. As already 
discussed, they are locally isomorphic to a tree, but they contain random large-scale loops which are not consistent 
with any crystalline structure. In particular they do not allow for a coherent partition of the original graph into two 
subgraphs as discussed in the context of the crystalline solution, the system is frustrated. So far, for these graphs we 
have only the liquid solution which, however, is in contradiction to numerical results for large /i. 

We therefore have to consider the possibility of the existence of a glassy high-density phase, or, in technical terms, 
we have to search for a replica symmetry broken solution (RSB). On the level of one-step RSB (IRSB) this can be 
achieved most easily using the cavity approach pj| . as developed for finite-connectivity systems in 12], and used also 
in for a similar lattice gas model. 

We first discuss a homogeneous IRSB solution, which is valid in particular for lattices not allowing for crystallization. 
At the end of this section, we also show shortly which modifications have to be made in order to account for the 
crystalline glass phase which was already mentioned at the end of the last section. 



A. The IRSB cavity solution 

We want to construct a non-periodic solution of Eq. 0. For small fj,, the lattice gas is very diluted, so the liquid 
solution is expected to be correct. For higher densities, the frustration in the systems becomes, however, relevant. In 
this regime, Eq. © is expected to have an exponential number of solutions, each one describing a thermodynamic 
state of the model. These states are, in particular, not spatially homogeneous. Let us therefore assume that the 
number A/'Ar(a;) of states of density uj of the grand-canonical potential satisfies 

AAwH-e^^("), (20) 

where S(w) > is the so-called complexity (or configurational entropy). It is assumed to be a growing concave function. 

We consider now a rooted tree with root jl as shown in Fig. O for j — 0, I = 0, and assume that the solutions 
for the cavity field x^' and the potential a;-'' are distributed according to some probability distribution W^{x^\bj^^). 
Repeating the iteration step represented in Fig. |21 we can calculate the distribution for the root from the distribution 
of the neighboring roots. 



i?"°(a;°",c^"°) = I \{\{dx^'du,''W\x^\u:'')5{x'''' ~x{x^\...,x''n) 
j=ii=i 



X S - ^ f E E + Af^itc.(x", ...,x'n)]- (21) 



^00 

^J=l 1=1 

We have used the number N^'' of the vertices in the corresponding rooted tree, i.e. 

fc p 



= E E + 1 • (22) 

3=1 1=1 



The change in the grand canonical potential is calculated via 

-fj,Afiitai partition function after iteration 



e 



partition function before iteration 

"00 I ^00 



n,tinLiN'+si') 



uuniii^i'+^ii 



k -I— r?) /- . ..jN ' (^^) 
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Note that Eq. (|21|l assumes factorization of the joint distribution R{x^^,uj^^; . . . ; x'''^, oj^p) of aU neighbors of vertex 
00. This would obviously be correct on a tree, where different rooted subtrees are disconnected and thus statistically 
independent. On a (generalized) random regular graph, these vertices are not the roots of kp disconnected rooted 
trees, so the factorization is not a priori true. Loops have, however, a typical length of ©(In TV), so once their short 
connection via vertex 00 is disregarded, they are far away. On the level of IRSB we therefore assume that the joint 
distribution factorizes, and we can use Eq. (|21|l . 

On these graphs, regular structures like a crystalline one are forbidden by the existence of random loops. We 
therefore assume homogeneity on the level of the field distributions, = R^^ = . . . = R^p = R, thus Eq. f?n 
becomes a self-consistent equation for R{x,lu). In order to solve this equation, we use assumption (|20|1 and write 

Rix, uj) ~ e^^('^) P('^) (x) , with / dxP^'^^ (2:) = 1 , (24) 



where p(") is the distribution of cavity fields with fixed co. In order to derive an equation for this function, we fix the 
density of the grand-canonical potential at some arbitrary value wo, and introduce the parameter 

to(wo) = - 3^(^o) (25) 

and linearize the complexity in a small vicinity V^^g of ujq: 

T,{lu) ~ S(wo) + rnii{uj — uq) . (26) 

Since S was assumed to be concave, there is a one-to-one correspondence between ujq and m. We define the 
distribution 

P(")(a;) = p^/ dujP^^\x). (27) 
Plugging this into Eq. H21|l . we find in the limit of a vanishing vicinity Vi^„ the self-consistent equation 

/k p 
Jl Jl da;J'F('")(a;^') 5{x - i(x", . . . , x'^^)) e-™MAJi.t„(xi\...,:.'='') (38) 
3=11=1 

which is called the factorized IRSB cavity equation for the cavity field distribution in our hard-sphere lattice gas 
model. 

To clarify the meaning of m, we calculate 

5(m) ^ ^ e"^"''"*"' = / dcjA6v(a;)e-^™^" - / dcj e^^^^")"™^'^) = g-JVm^<E>(m) (29) 

with the sum running over all states a, i.e. over all solutions of Eq. Note that, for m — 1, this reproduces the 
grand-canonical partition function. In the thermodynamic limit, <I>(m) is given by the saddle-point approximation, 

$(m)=w ^S(cj), with 9^S(cj)=m^, (30) 

m/i 

which reproduces the definition of m(cj) in Eq. ((23) • From Eq. ^ follows that 

d„Mm) = draiU + ^S(W) - —d^^{uj)draLO = ^S(a;) . (31) 

Let us compare Eq. H30|) with the usual equation for the entropy 

$(m)=w ^S(u;) u} = -p--s. (32) 

m/i /i 

There is an obvious similarity: As s counts the number of microscopic configurations, S counts the number of states. 
As /i serves as a control parameter selecting configurations of a given density p, the parameter m allows to focus 
on states of arbitrarily given lo. Therefore $(m) plays, on the level of states, the same role as the grand-canonical 
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potential on the level of configurations. Following this analogy, we can calculate $(to) as done in Sec. lIVBl for the 
grand-canonical potential: 

Hrri) = _^i^±i) A$iink(m) + A$,itc(™) (33) 
p+1 

with 

, k+l p 



and 



e 



j=i 1=1 



1=1 



Aflsite and AfJnnk are the direct generalizations of the corresponding quantities in Sec. IIVBI 



nL^'(i+^«)' n-;!nLi(i+^^o 



(36) 



1. Static and dynamic glass transition 

Let US consider the IRSB equation H28|) first for m = 1. As already mentioned, S(m — 1) then gives the usual 
partition function, 

5(/x) - J dLje^(^(")-'^") , (37) 

where the integration runs over the interval (wmin, Wmax) with positive S(a;). 

Let us first discuss the case p > 1. For small /i, Eq. (|28|l has only the liquid solution P{x) — 6{x — x*). The particle 
density is low enough that the frustration resulting from long loops is not relevant. At some chemical potential fi — fid, 
the first non-trivial solution appears discontinuously. This point is called the dynamical glass transition. For fi > fid 
the system has exponentially many states, the complexity S(a>s) > is positive. The grand-canonical potential ujs of 
one state can be calculated via the saddle point method from Eq. (|37|l . i.e. as a solution of /i = duj'S{uj). The potential 
uj* of the liquid solution is given by lu* — <i>(m = 1). With Eq. we conclude uJs > i^* ■ The states born at fid 
are therefore metastable ones, the thermodynamically stable solution remains still the liquid one. The dynamics of 
the system becomes, however, trapped already at fid by the exponential number of metastable states, i.e. the system 
remains dynamically out of equilibrium, cf. the discussion in 15]. 

Increasing fi beyond fid, the complexity Tt(ujs) decreases, until it vanishes at /x = fig. According to Eq. H3U|) . the 
grand-canonical potential lOs of the IRSB solution equals here the one of the liquid solution, and fi = fis indicates 
thus the location of a static, i.e. thermodynamic glass transition. 

For fi > fis the naive application of the saddle-point method leads to an unphysical solution a;^ < o^niin of negative 
complexity. The correct exponentially dominating contribution in Eq. (|37|l is thus given at the interval limit lUs = i^min- 
In this situation we need a "corrected" saddle-point equation m/i — 9ijl](u;) (cf. Eq. H3U|) '). which defines the value 
rris = -duj^iuj = LUs). According to Eq. 13UI) we thus obtain the grand-canonical potential of the IRSB equilibrium 
via ^(nis) = LUs = "-^min- Using S(ciJs) — ^(cjmin) = and Ea. (|31|) we thus find 

a„,$(TO = m,) = 0. (38) 

Also other values of m have a physical significance. For m < we get metastable states with lu > lOs- These are 
important for understanding the dynamics of the system. A particular importance has here the value LUd where the 
complexity is maximized, cf. Sec lVIl Larger values m > Wg lead to negative complexity which can be reinterpreted in 
the context of atypical graphs - their complexity corresponds to the probability that a (generalized) random regular 
graph has states of this smaller grand-canonical potential. An example for such an untypical graph is the multi-partite 
one allowing for a crystalline solution as discussed before. 

The situation is easier for p — 1, i.e. for standard Bethe lattices. Here fid and fis coincide and mark a continuous 
spin-glass transition. The region, where the dynamics is dominated by an exponential number of metastable states, 
disappears. 
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2. The spin-glass instability of the liquid solution 



An upper bound for /i^ is given by the spin-glass instability of the liquid solution. This instability means that a 
small perturbation of x* is amplified during the iteration of Eq. H28|l . The criterion for the instability is 



kp 



dx{x* , . . . ,x*) 



dx^ 



> 1. 



(39) 



Equality holds at critical values of the chemical potential, and pg of the density. According to Eq. © this 
happens at 



1 



Pa 



/kp+1 



(40) 



Above this point, the liquid solution is locally unstable. 
For p = I, also this instability coincides with fi^ and fi^ 



3. Solving the IRSB equation via population dynamics 

Solving Eq. (|28|l seems to be too complicated, we therefore have to solve it numerically. The main idea hereby is 
to iterate the self-consistency equation until a fixed point is reached. This happens by representing first p(™) by a 
finite but large population {xi, . . . ,xm}, M ^ 1 of cavity fields, and by updating this population via a population 
dynamical algorithm In every step kp fields are selected randomly from the population, and a new field is 

calculated by x. This field replaces one or more of the old fields in the population, where the number of substituted 
fields is a function of the reweighting factor in Eq. H28|l . 

For our example A; = 3, p = 2, we find the following results: From Eq. (|40|) we conclude that pg = — 0.290. It 

is, in particular, larger than pkrit — -j^p^ — 0.25, which is true for all k > p. With Eg. I13l we find pg ~ 5.89, this value 
is confirmed by the population-dynamical algorithm. 

Further on, we find that at pd — 4.85 the first non-trivial solution for p(™) emerges discontinuously. In principle 
also the static spin-glass transition ps can be determined by the population-dynamical algorithm. The problem is, 
however, that it results from a comparison of the grand-canonical potentials of the liquid and the glassy solutions, 
which are subject to large fluctuations. Even a large statistics allows only for a rough estimate of ps ~ 5.0. 

Note that, for the case considered here, the dynamical and the static RSB transition are located beyond the 
inverse-crystallization point, where the liquid solution becomes locally unstable. It is, however, possible to validate 
the correctness of the result by numerical simulations which have to be performed on a (generalized) random regular 
graph. These graphs have a large-scale loop structure which is not consistent with any crystalline packing, cf. the 
discussion in Sec.^ and thus the system undergoes directly the glass transition if compactified until pd- As can be 
seen in Fig. O the compaction dynamics falls out of equilibrium even before the dynamical glass transition, but slower 
and slower compaction allows for a closer approach to pd, showing thus the compaction-rate effects which are typical 
for glass formers. 



4- Glass-instability of the inverse crystalline phase (k = 3, p — 2) 



A stability analysis of the inverse crystalline phase toward a IRSB "glassy" phase requires to embed the inverse 
crystalline solution in the IRSB formalism. To this purpose we rewrite the IRSB equation H28f) without assuming 
homogeneity of the field distributions P"", P^^, . . . , P'^P; 

P°"{x) (xT[P^\...,P''P]{x) (41) 

where contains the r.h.s. of the IRSB equation. The dependence on m is eliminated by putting m = 1 which can 
be justified as we are only interested in transitions points to the glassy phase. The canonical embedding now reads 

P(°) cx .F[P(1) , P(°) , . . . , P("' ;...], P(i) oc .F[P(°) , . . . , P(°) ;...], (42) 

where the displayed blocks in the argument of J-' must be repeated k times similar to Eq. ||^. The distributions 
P(0) and P(i) are simply (5-functions at the values x^^'' and x'^^^ of the local fields according to the inverse crystalline 
solution. The local instability of the iteration defined by Eq. can be detected by evaluating for P^"^ and P^^^ 
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FIG. 8: Compactification experiments for a generalized random regular graph fc = 3, p = 2, and A'' = 450000 sites. The full 
lines give the analytical results for the liquid and crystalline densities, whereas the symbols correspond to compaction rates 
SfJi/5t = 0.1/10 MCs (circles), 0.1/100 MCs (squares) and 0.1/1000 MCs (diamonds). Obviously, the systems falls out of the 
liquid equilibrium when the equilibration time starts to exceed the waiting time at the corresponding chemical potential, i.e. 
before the latter reaches jid- 



chosen as Gaussian distributions with respective variances and e'"'^'. Under neglecting of the reweighting factors 
in T one easily calculates the variances e*^"^ and e*^^^ after one step of the iteration: 



aT(a;(i),a;(o),...,a;(o);...) 




5i(a;(i),a;(0),...,a;(0);...) 



AO) 



(43) 



We note that the equations dehver precisely the spin-glass instabihty of the liquid solution when we choose all local 
fields equal to x*. In the case of local fields a;^"^ and a;^^^ a stability criterion can be derived by evaluating the 
eigenvalues of the Jacobian matrix attributed to the iteration defined by Eq. (|43|) . We find that for the inverse 
crystalline solution e^*^^ — e*^^) = is a stable fixed point of the iteration only if /i < Hg^inv — 5.31 which is smaller 
than ~ 5.89. We can check this result via population dynamics where we encode the distributions P^^^ and P'-'^' 
using two populations of fields. It turns out that the local instability Hg^mv is exact and furthermore we can show that 
for /Lt < fj.g^inv only the trivial field distributions corresponding to the crystalline and inverse crystalline solutions solve 
the IRSB equations 1421) . This is obviously in contradiction with the results from section lTV El where the existence of 
a glassy crystalline phase is shown for fi > /icg ~ 4.1. We can thus conclude that the IRSB formalism is not sufficient 
to describe the onset of a glassy phase. 

A way out of this dilemma might consist in passing to a description in 2RSB which would mean to work with two 
distributions of field distributions attributed to the sub-lattices. A 2RSB treatment of the problem seems, however, 
complicated to handle. For the moment, the construction of a "glassy" solution which would illuminate the discrepancy 
between ^g^mv and ^cg remains an open problem. 



VI. CLOSE-PACKING LIMIT (^i ^ oo) 

The IRSB solution developed in the preceding section FV Al cannot be calculated analytically in the case of general 
jj, where a solution is obtained by a numerical method (population dynamics). 

For the limiting case /i — s- cx) which corresponds to a close-packing, however, one can easily show that the transformed 
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local fields 



1 ■ r..' - 1 



= - \n{x^' + 1) = _ In ' * (44) 



must have a particularly simple distribution in the IRSB formalism which enables us to undertake an explicit calcu- 
lation. 

We first note that for fi oo the grand-canonical partition function is dominated by the terms that correspond to 
configurations with a maximal number A^^y^ of particles, i.e. S^y^ oc e^^"/' asymptotically with a factor that counts 
the number of close-packings and docs not depend on fi. This in mind, we find that Eq. 1441 becomes 

h^^ = max(iVj' , Ni^) - TVf (45) 

which implies that the local fields h^^ must be nonnegative integer numbers. 

Using again the properties of the partition function for n oo we can also rewrite the iteration equations 01 in 
terms of maximal numbers of particles: 



N 



00 



^max( f]7Vi', [Ni' -I- iVi™ | / = 1, . . . ,p} J (46) 



k p 

= 1 + Xl£^e'- (47) 

j=l 1=1 

We use these relations to determine the maximal value of putting j = and / = without loss of generality. 
Following 145|l the interesting case is N^'^ > N^^ where we find 



hoo ^ nOO _ ^00 



p 



k p 



\ 



' Ni 



J 



= l \l=l m=l / i=l i = l 

k 

= 1 - E max(0, {Nl' - iV^' | / = 1, . . . , p}) 

i=i 

which implies h^^ < 1. Thus we have shown that the transformed local fields /i-'' satisfy h^^ £ {0, 1}. 

This result, which shows the interest of working with /i-'' instead of x-'', enables us to reduce drastically the degrees 
of freedom in the IRSB field distribution by making a simple one-parameter ansatz: 

Pih)^Pod{h)+piS{h-l) (48) 

where normalization requires Po + Pi = 1. 

The pair of parameters po,Pi must be chosen such that the ansatz solves the IRSB cavity equation (|28|l which can 
be directly rewritten for the transformed local fields: 

F('")(/iO") cx f ]^f[d/i^''F(")(/i^'')(5(/i°"-/i(/ii\...,/i'=P))e-"'''^^--(''"--''^^^ (49) 

j=ii=i 

where the recursion relation h of the local fields and the change in the grand-canonical potential ASlitor per iteration 
step are now given by 

^^Hh^\...,h^n = 1+ (50) 



and 



nL(i-p+Ef=i 



e 



p ^fih' 



(51) 
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as follows directly from Eq. © and Eq. I|23|l . 

For practical reasons we define a vector V with k components that encodes the relevant information of a local field 
configuration {h^^, . . . , h'^P}: 



i^, = Y.h='. (52) 

1=1 

As h^^ satisfies h^^ E {0, 1}, i^j counts the non-vanishing fields in the jth clique. Furthermore we define a norm by 

\iy\=#Wj^O\j = l,...,k}. (53) 

We can see from Eq. (fSU)! and Eq. that h and Afijter are functions of only. We calculate these quantities 
explicitly: 



V n-=i(i-p + ELie^'^^'), 



and 



1 if|i7|=0 
if \V\ > 



Ariit^- = - In ' 



(54) 



^ V^■=l(l-P + ELle'^'^^')+e^ 

k p 



3 = 1 1=1 ' \j=l 1 = 1 

k if'' 

^i/,--ln Hil - + J^je'^) + 

3 = 1 ^ \3 = 1 

-1 if \v\ = 

= E|=i max(0, - 1) 
= niax(l, Vj) - fc if > 



(55) 



These results permit us to evaluate the IRSB Eq. 149|) for the ansatz (|48f) as they reduce the integrations in the 
equation to a simple combinatoric problem. 

We see from Eq. (|54f) that only the configuration v — Q contributes to /i™ = 1. With the ansatz for P{h) we deduce 
immediately from the IRSB equation that 

Pi (X (56) 



holds, where we have set y = m/i and used AfJjtor from Eq. (|55|l . 

In a second step we must sum up the contributions to hP'^ = that come from all configurations with v ^ Q. To 
this purpose we note that there are 



k 

P 



n 

3 = 1 

different realizations of a configuration with fields {h^^, . . . , h'^P}. 
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For each configuration we collect the probabilities emerging from the ansatz yielding a factor 

Pi Po 

In addition, we have to take into account the reweighting factor 

g-yAfiitor _ gj/fc g-yl2 j=i max(l,i/j) ^ 

We obtain the contribution from a configuration v hy taking the product of the three terms 

k 



n 

i=i 



The overall contribution to /i*^" — can now be calculated by summation over the components of i7 ^ 0: 



ui=0 i/fc=0 j = l 



correction for 



j=l \ uj=0 ^ J 



{ey{p^e-y + po)" -e>g + pg} " pI" ■ (57) 



correction for 

Making use of normalization (po +Pi = 1) we end up with the iterative equations for the probabilities po and pi: 

. , . {ey{pie-y+P,r-Pl{ey^\)t-p\^ 

Pq{PO:Pi) = -T 1 (58i) 

{evipie-v + po)P - p^oi(^y - 1)} +Pa''{ey - 1) 

kp y 

Pi{Po,Pi) = ^ -j: 1 . (58ii) 

{evipie-y + po)P - p^oi^^y - 1)} + Policy - 1) 

The equations can be solved numerically for given y with the help of a computer algebra system. 

It remains the question which value y = ys we must choose in order to describe thermodynamic equilibrium. The 

answer is given in section FVAI via the condition that the configurational entropy must vanish, i.e. S(?/s) = 0. We use 
Eq. I^with the substitution y = m/K to find 

E(y) = y'dy<^>{y) . (59) 

which provides = ?/s) = as an alternative condition for y^. 

Consequently, we shall compute $(y) to complete the picture of the close-packing limit. $(y) is given by Eq. H33() 
with Eq. H34|) and H34|l as well as Eq. (|36|l . We proceed analogously to the derivation of the equations for pQ^pi: in a 
first step, we evaluate Eq. H36|l with transformed local fields in the limit ^ — > oo. In a second step, the integrations 
in Eq. H34|) and H34|) are performed. This reads: 

P+i -, /p+i \ 

1=1 ^ \i=i 



link = > , n^ ' in > e'-- - p 



V - iln(l - J/ + ve^) 



[I — >oo 



V — \ 11 V > \) 
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where we have set v = X]f=i ^^^^ ■ Further on we have 

fe+i 



' fc+i 



j=i 1=1 ^ \ 1=1 
fe+i ^ /fc+i 



if = 



J2';tl max(l, i/j) - (fc + 1) if W\ > 



(61) 



Note that the expression for Arigitc is obtained from Afiitcr by the substitution k k+1. We perform the integrations 
to obtain 



and e-^^*='- = {e^Pie-'' + Pof - pg(e^ - 1)}''^' + P^'''^'^''(e^ - 1) , 

where the latter term also resuhs from the sum of the contributions (|56|l and (|57|l after substituting k ^ k + 1. 
The knowledge of $(?/) can serve to access the particle density p: 

p = d^\nE = -d^ifi^) = -dyiy^). 

SpeciaUy in thermodynamic equihbrium the particle density is given by 

as dy^ vanishes at ys- 



(62) 
(63) 

(64) 
(65) 
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P 

FIG. 9: IRSB solution for the complexity E as a function of the particle density p in the limit /i — » oo (fc = 3,p = 2). S 
being a concave function, only the full line has physical meaning. E < is attributed to untypical graphs that have vanishing 
statistical weight. 
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Eq. H64I) together with Eq. (|59|) deliver the complexity as a curve parameterized by y (Fig. |2J|. The complexity 
S has a maximum at pd < Ps while it vanishes for the equilibrium density ps- The meaning of pd for the dynamics of 
the system, in particular for dynamics induced by local algorithms, is discussed in the subsequent section. 

Some results for the particle density of a close-packing poc = Ps and the density pd with maximal complexity as 
obtained in the IRSB solution are shown for different k and p in Fig. 1101 

A. Comparison with numerical bounds 

We have compared the analytical IRSB results for p^o and pd with numerical lower bounds for p^o that were 
obtained in Monte-Carlo simulations. For given k and p random graphs of size N = {p+1) ■ 500 000 that have locally 
the structure of a generalized Bethe lattice were generated. Roughly speaking, this can be achieved by iterating 
the following procedure until every vertex has (fc -f l)p neighbors: choose randomly p + 1 vertices with less then 
(fc -|- l)p neighbors and fully connect them to form a clique. Furthermore we introduce local dynamics on the lattice 
by applying the following local rules to a randomly chosen site: (i) if the site is empty: introduce a particle, (ii) if 
the site is occupied: move the particle to a randomly chosen neighboring site. If the application of a rule conflicts 
with the hard-core constraint, the rule is not applied. The time 1 MCs corresponds again to choosing N times a site 
at random. Fig. 1101 shows the particle density that the initially empty system {t = OMCs) has reached at MC-time 
t = 100 000 MCs. The error can be estimated to less than 0.1 % by comparing the results for different realizations of the 
graph. The asymptotic behavior of the density as a function of time can be characterized by the following observations. 
The density at t = 50 000 MCs {t = 75 000 MCs) is about 99.96% (99.99%) of the density at t = 100 000 MCs for 
small p and about 99.86 % (99.95 %) for great p. These figures illustrate that the compaction process becomes very 
slow in the second half of the simulated time frame. We do not expect the system to equilibrate on a finite time scale 
due the presence of glassy states. The quality of the numerical results, which constitute lower bounds for the exact 
values of poo , is thus determined by the finiteness of the simulated systems and the simulated times which were chosen 
maximal within the limits of reasonable computation time. 




1 2 3 4 1 2 3 4 

p p 

FIG. 10: Left: particle density pca of a close-packing (-I-) and density pd with maximal complexity (x) from the IRSB solution 
compared with numerical lower bounds (0) for poo (from top to bottom: k — 1,2, 3, 4). For k = 1, poo and pd coincide. Right: 
the difference of pd (x) and the numerical bound (0) to poo = Ps is shown for k = 2,3,4 (top to bottom). 

We notice that the IRSB results for poo are never violated by the numerical bounds. In the case fc = 1 the agreement 
between numerical and analytic results is very good which can be understood e.g. assuming p = 1 where the "random" 
graph forms simply a circle which can be always packed in the crystalline ground state up to a misfit if N is odd. For 
k > 1 the numerical bound stays always below poo with a gap that grows with increasing p up to 4.0 %. On the other 
hand, we find a good agreement of the numerical bounds with the results for pd. We shall explain this observation 
with a remark on the configuration space structure as obtained within the IRSB solution. 
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We have assumed as a starting point for the derivation of the IRSB cavity equation that a non- vanishing complexity 
is equivalent to the presence of an exponential (in N) number of metastable states in the configuration space {ergodicity 
breaking). As different metastable states are only connected via configurations of lower densities any local algorithm, 
trying to increase the number of particles in every step, gets trapped in a metastable state. Using the definition of the 
complexity (|20|l it becomes clear that the state in question will have almost always a particle density p ~ < poo 
which maximizes the complexity (cf. Fig. |5J. The numerical bounds that result from a local algorithm give indeed 
strong indication that the concepts of the IRSB formalism have a very concrete meaning for algorithmically induced 
as well as physical dynamics. 

Note that there is recent work by Montanari and Ricci-Tersenghi who performed numerical cooling experiments for 
p-spin glasses (16} , a model which is somewhat similar to our lattice gas model. The behavior for different cooling-rates 
underlines the relevance of the so called iso-complexity curves. In this context, the numerical experiments performed 
in our study constitute a limiting case of very high compaction-rate which would lead to densities significantly below 
Poo if the particle movement, which has no equivalent in the spin-glass context, was omitted. An investigation for 
different compaction-rates might give interesting complementary insights concerning the glassy phase of the lattice 
gas. 

We conclude with the remark that there is a priori no indication for the correctness of the results that were calculated 
under the restrictions of the IRSB ansatz. In the next section we shall sketch the baselines of a stability theory for 
IRSB in order to check the reliability of the preceding results. 



B. Stability analysis of the IRSB solution 

In this section we give a brief summary of concepts concerning the local stability analysis of a IRSB solution towards 
a 2RSB solution. A detailed discussion of the topic was given by Rivoire et al. J3] for the case a Bethe lattice gas 
which shares common properties with the system under study. 

We start with the introduction of the 2RSB formalism. Compared to the IRSB solution which was based on the 
organization of configurations in states, 2RSB acts on the level of states by assembling states in clusters to which we 
attribute grand-canonical potentials lOc- The number of states within a cluster is given by referring to uJc'- 

NsiuJs) ^ e''y^^^'-^^\ (66) 

where ys must be understood in analogy to Eq. H25|l as the slope of the complexity of the states within the cluster. 
Using this relation as a starting point, the IRSB cavity equation can be derived with the same arguments as in section 
IV Al yielding 

/k p 
j=i 1=1 

when we discard the homogeneity assumption of the distributions P^K The shape of the equation invites to iterate 
the transfer which starts with the recursion x{{x^'^}) to deliver the IRSB equation through the introduction of field 
distributions. Concretely, this means to start with the recursion P[{P-''}] and assume distributions of field distributions 
(clustering). As in the IRSB derivation we introduce the number of clusters with grand-canonical potential uJc with 
respect to a reference point ujq: 

J\fc{tOc) - e^''^ (68) 

The clustering is realized by the introduction of distributions Q[P] for the field distributions. Using the arguments 
from section FVAI the so-called 2RSB cavity equations can be derived: 

/k p 
Jl J|l?P^'Q[pj']<5(P-P[{P-''}])e-''=^*""[{^">l . (69) 
j=i 1=1 

where A$itcr is the change in grand-canonical potential corresponding to one step of the iteration P[{P^'}]. It can 
be obtained by evaluating AOitcr with respect to the distributions P^': 

g-y.A*,„[{P^'}] ^ f f[f[dx^lpjl(^^ll)e-y'^^^^-^{-''^l (70) 
•' 3 = 11=1 
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The hierarchical scheme of distributions can be continued to higher levels of RSB, leading to full replica symmetry 
breaking (frsb) in the limit of infinite continuation. A proof for the correctness of a IRSB solution would require to 
show that every description in higher levels of RSB coincides with the IRSB description. This is of course not feasible 
if only the IRSB solution is explicitly known. In this case one can show the local stability of the IRSB solution 
toward 2RSB which is believed to be a good indication for the correctness of the IRSB solution as no example for a 
discontinuous transition from IRSB to RSB of higher levels is known so far. 

It was shown be Montanari and Ricci-Tersenghi that there are two possible kinds of instabilities that may occur |19| | . 
These result from two different descriptions of a IRSB solution within the 2RSB formalism. One way is to assume 
the existence of a single trivial cluster where the states are given by P*{x), with P* being the IRSB solution. Thus 
we have Q[P] = S[P — P*] and the corresponding instability can be detected by making the ansatz Q[P] = f[P — P*] 
where / is a functional with support around the null function. This corresponds to an assembling of the states that 
organize in clusters. The process is called aggregation of states and the instability is classified type I. The second 
possibility consists in encoding the field distribution P* (x) in the cluster structure while the states are trivial. For 
discrete P* we write P* = J2aP<>-^<^ '^i^h 6a{-) = 5{-—a). The cluster structure is then given by Q[P] — X]aPa<5[P— ^a] 
and the instability is detected by Q[P] = ^aPafa[P — So] with fa having support around the null function. At the 
instability each state can be considered as a germ giving birth to a cluster. The process is called fragmentation of 
states and the instability is classified type II. 

We shall not stress on the details of the concrete calculations which can be found in the literature for various models 
including combinatorial problems ITtI 111 IT9ll20ll2lll2^ . 

In order to investigate the stability of the IRSB solution in the close-packing limit — > oo) we focus on the type-I 
instability. We note that there is a strong analogy between the calculation of the spin-glass instability, which was 
detected by putting P{x) = f{x — x*) where / is a smeared (5-function, and the ansatz Q[P] = f[P — P*] suitable to 
find the type I instability. By transporting the calculations for the spin-glass instability to the next higher level in 
the RSB hierarchy the stability criterion fcp|Amax(?/)P < 1, which reads precisely like Eq. (|39fl . has been derived by 
Rivoire et al. Qj. Here Amax(y) is the eigenvalue of greatest modulus of the matrix 



which is the Jacobian matrix associated with the iteration of the IRSB field distributions Eq. (|58|l . 

A computation of Aniax(y) yields that we have kp |Ainax(2/)P > 1 at y = as well as y — for all the cases shown in 
Fig. ^1 i.e. the IRSB solution is locally unstable toward 2RSB. We can conclude that the results for pd and poo that 
were derived in the IRSB context are not exact. Following the general belief that the instability of a IRSB solution 
indicates frsh, this would even mean ps — pd for /i — > oo, which is not suggested by the numerical findings. It must 
however be considered that a system cannot be brought instantaneously to infinite /i via numerical compaction. Under 
the action of a MC-algorithm the simulated system will pass through a chain of genuine non-equilibrium states. By 
mapping these on equilibrium states according to their density we can approximately attribute a chemical potential 
li{t) to the system. It might seem more adequate to rely on a stability analysis of the IRSB solution for = /x^j in 
order to decide if the compacted system is governed by metastable states. Unfortunately, this analysis seems to be 
out of reach for present analytical techniques. 



To conclude, we have investigated in great detail the thermodynamic behavior of a hard-sphere lattice gas model 
on generalized Bethe lattices. 

On usual Bethe lattices, which do not contain (short) loops, the model shows a very simple behavior: At low density, 
the system is found in a liquid phase, and the local density is homogeneous. At a certain point, the system undergoes 
a second-order freezing transition. The character of this transition depends on the global graph structure: The frozen 
phase is either crystalline (if long loops of odd length are almost absent in the graph), or it is given by a spin-glass 
phase (if the Bethe lattice is defined as a regular random graph, and the crystalline structure is inconsistent with the 
large-scale structure of the graph). 

The introduction of short loops changes this behavior drastically. We have therefore studied generalized Bethe 
lattices, which are locally characterized by the existence of many short loops, but on a coarse-grained level form 
locally tree-like hyper-graphs. This mixed structure allows for an analytical solution of the hard-sphere lattice gas. 

Again, the low-density phase describes a spatially homogeneous liquid. If compactified, the system undergoes either 
a first-order crystallization transition or a discontinuous glass transition which is characterized by a broken replica 
symmetry, and by the existence of an exponentially large number of metastable states. We also find more exotic 
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high-density phases which are at least metastable and can thus be found also in numerical simulations. The first of 
these phases is an inverse-crystalline one. Its sublattice densities are inverted, as compared to the densest crystalline 
packing. For even higher density, this phase becomes unstable with respect to a crystalline glass transition: Even if 
the local densities are frozen to random values, their distributions are organized in a regular, i.e. inverse crystalline 
manner. 

This means that hard-sphere lattice gases on generalized Bethe lattices may serve as valid, microscopically mo- 
tivated and analytically tractable mean-field models for the glass transition. Do they, however, also present good 
approximations to the behavior in low dimensions? The hard-sphere model is in fact well studied on various two- and 
three-dimensional lattices, the most famous cases are the hard-hexagon model and its solution by the corner-transfer 
matrix |23.. 24, 25], and the three-dimensional cubic lattice (26i | which can be analyzed by means of series expansions. 
These two cases show second-order crystallization transitions. In we have, however, introduced a finite-dimensional 
lattice mimicking the structure of the generalized Bethe lattice, and this model shows in fact a glass-like slowing down. 
The latter is, however, much less pronounced due to the existence of activated processes in finite dimensions. 

At this point, another remark concerning the relation of the generalized Bethe-lattice models to p-spin glass models 
is necessary. The latter model class was extensively used in particular in the last decade in the context of the 
glass transition |^|^. On one hand, they allow for an analytical treatment of the model statics and dynamics, the 
latter being formally equivalent to the mode-coupling approach to more realistic glass-forming liquids, p-spin glasses 
are, however, completetly missing a microscopic motivation. One point here is the existence of unphysical multi- 
spin interactions. If we look, however, to our generalized Bethe-lattice model, the introduction of cliques effectively 
introduces multi-site interactions - only one site per clique can be occupied - even if these are defined via groups 
of pair interactions. This allows to reinterprete also the p-spin interactions in spin-glass models as a kind of coarse 
grained effective interaction. 

All the considerations in this paper are restricted to the thermodynamic, i.e. equilibrium behavior of the model. In 
particular in the glassy regions, we know that the dynamics becomes very slow, and the system is practically always 
out of equilibrium. The analytical description of the dynamics of finite-connectivity models is, however, a much more 
involved problem, and it is to a large extent an unsolved one. In a subsequent publication, we will confront our findings 
on the static behavior with an approximate analysis of the dynamics [l^ . The approach taken there is related to a 
projective approximation scheme introduced in 27, 28] in the context of algorithms for combinatorial optimization, 
and then applied to simple Ising ferromagnets in j29j . Within this scheme, the dynamics is projected to the dynamics 
of a small set of global observables, and the resulting equations are closed on the basis of a pseudo-equilibrium ansatz 
in an enlarged ensemble. The selection of this generalized ensemble, and the technical realization of the projection 
are based on the tools developed in this article, and therefore will form a natural continuation of the present work. 
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